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The velocity fluctuations present in macroscopically homogeneous suspensions of neu- 
trally buoyant, non-Brownian spheres undergoing simple shear flow, and their depen- 
dence on the microstructure developed by the suspensions, are investigated in the limit 
of vanishingiy small Reynolds numbers using Stokesian dynamics simulations. We show 
that, in the dilute limit, the standard deviation of the velocity fluctuations (the so-called 
suspension temperature) is proportional to the volume fraction, in both the transverse 
and the flow directions, and that a theoretical prediction, which considers only for the hy- 
drodynamic interactions between isolated pairs of spheres, is in good agreement with the 
numerical results at low concentrations. Furthermore, we show that the whole velocity 
autocorrelation function can be predicted, in the dilute limit, based purely in two-particle 
encounters. We also simulate the velocity fluctuations that would result from a random 
hard-sphere distribution of spheres in simple shear flow, and thereby investigate the ef- 
fects of the microstructure on the velocity fluctuations. Analogous results are discussed 
for the fluctuations in the angular velocity of the suspended spheres. In addition, we 
present the probability density functions for all the linear and angular velocity compo- 
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nents, and for three different concentrations, showing a transition from a Gaussian to 

an Exponential and finaUy to a Stretched Exponential functional form as the volume 

fraction is decreased. 

The simulations include a non-hydrodynamic repulsive force between the spheres which, 
although extremely short ranged, leads to the development of fore-aft asymmetric distri- 
butions for large enough volume fractions, if the range of that force is kept unchanged. 
On the other hand, we show that, although the pair distribution function recovers its 
fore-aft symmetry in dilute suspensions, it remains anisotropic and that this anisotropy 
can be accurately described by assuming the complete absence of any permanent doublets 
of spheres. 

We also present a simple correction to the analysis of laser-Doppler velocimetry mea- 
surements, that only takes into account the mean angular rotation of the spheres in the 
vorticity direction, and which substantially improves the interpretation of these measure- 
ments at low volume fractions. 



1. Introduction 

The problem of determining the velocity fluctuations in suspensions of non-Brownian 
solid spheres in Stokes flows is one of long-standing difficulty due to the underlying long- 
range many-body hydrodynamic interactions between the suspended particles. Even an 
apparently very simple case, that of determining the dependence on the shear rate of 
the velocity fluctuat ions in simple shear flows, remains a matter of some controversy 



i Shaplev et al 



2002f) . What is clear is that, although the suspension might be homoge- 
neous at macroscopic scales, the continuous rearrangements in the suspension microstruc- 
ture and the corresponding hydrodynamic interactions between particles lead to fluctu- 
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ations in the particle velocities about their mean values, in both the transverse and the 

flow directions. 



I 

In our previous work ^Drazer et al. 



2OO2I to be referred hereafter as paper I) , we showed 
that the dynamics of sheared suspensions is chaotic and offered evidence that the chaotic 
motion is responsible for the loss of memory in the evolution of the system. This loss of 
memory, coupled with the fluctuations in the velocity of the spheres, ultimately leads to 
the phenomenon of shear-induced particle diffusion. 

The variance, or the standard deviation (STD), of the velocity fluctuations is the sim- 
plest measure of the magnitude of such fluctuations and is sometimes referred to as the 
suspension temperature, which in the case of an anisotropic motion of the suspended 
spheres would actually be a tensor (covariance matrix). The suspension temperature is 
relevant to the migration and diffusion of particles in shear flows, phenomena that occur 
in a wide variety of natural as we ll as engineering pr oblems, ranging from t he dispersion 



and migration of red b lood cells IjBishop et al 



2(M 



Cotz et al 



20021) to the food industry l|Cullen et 



2flfl3r . hence it is important to determine its properties. In particular, 
we are interested in the dependence of the velocity fluctuations on the concentration 
and microstructure of the suspension. Unfortmiately, and in contrast to the well-studied 
sedimentation problem, velocity fluctuatio ns in sheared su s pensions have receiv e d littl e 
attention thus far. In re cent experiments. 



Lvon fc Lea 



{12m and 



Shaplev et al 



Averbakh et al 



Shaulv et al 



(1992D; 



( 2002() used laser-Doppler velocimetry (LDV) to 



measure the vel o city fl uctu ations in conc e ntrate d suspensions of monodisperse spheres. 



Averbakh et al 



{IME) and 



Shaulv et al. 



1 199711 measured such velocity fluctuations in 



rectangular ducts and found that the STD's, both along and transverse to the flow, 
depend linearly on the shear rate (or on the maximum velocity inside the rectangular 
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channel), as expected in the Stokes hmit f 



'hubiiF wild A. 



Lvon fc Lea: 



■Act 



vos 



1 1998r measured the time- 



averaged local STD in the direction of the flow, for concentrated suspensions flowing 
in a two-dimensional rectangular channel , and also found a linear dependence on the 



volumetric flow rate. 



Shaplev et al 



(2003) presented the first detailed measurements of 
the velocity fluctuations in both the transverse and the flow directions, as well as of the 
dependence of the suspension temperature on the volume fraction and shear rate, for 
suspensions undergoing simple shear flow in a Couette device. Their results stress the 
difficulties encountered in such measurements and the discrepancy among different exper- 
imental results. They found a highly anisotropic temperature tensor, with the magnitude 
of the fluctuations in both transverse componen ts of the velocity smaller than that in the 



direction of the bulk flow. 



Shaplev et al. 



(2003) ^Iso found that the temperature is not 



monotonically increasing with volume fraction, as is usually expected, but shows a differ- 
ent behavior for each of its components. Specifically, the component of the temperature 
in the direction of flow was found to decrease with concentration, that in the direction 
of the gradient stayed constant, while that in the vorticity direction initially increased in 
magnitude with increasing concentration s and then dec r eased for concentrations larger 
than 40%. Finally, and most surprisingly. 



Shaplev et al 



1 2002^ found that, whereas the 



fluctuations in the direction of the flow increased linearly with shear rate (as expected for 
any flow in the Stokes regime), the STD in the vorticity direction increased non- linearly, 
while that in the gradient direction slightly decreased with shear rate. 

In terms of the suspension spatial structure, although a larger body of experimen- 
tal information exists concerning the microscopic structure developed by suspensions of 

I Let us note that in these experiments, the main contribution to the measured STD's were 

not actual fluctuations but migration and angular velocities of the suspended particles, as noted 
by the authors. 



Micro structure and velocity fluctuations 5 
monodisperse, non-Brownian spheres undergoing linear shear flow, no measurements of 

how the velocity fluctuations are affected by the suspension micr ostructure appear to have 



been conducted thus far. Recall that the experimental work of iGadala-Maria fc Acrivos 



( 198C1|) provided, for the first time, clear evidence that concentrated suspensions of monodis- 
perse, non-Brownian spheres develop an anisotropic structure when sheared. They showed 
that, when the direction of shear was reversed, the shear stress measured in a parallel 
plate device underwent a transient response not present when the shearing was started 
again in the same direction, and thereby concluded that the underlying structure was 
not only anisotropic but asymmetric under reversal of the flow direction, i.e. fore-aft 
asymmetric. Their oscillatory experiments showed similar results, in that the measured 
dynamic viscosity /i', although independent of the frequency of oscillation at low fre- 
quency, was consistently smaller that the shear viscosity /i of the suspension, stressing 
again t he presence of a microscopic structure induced by the shear. In recent experi- 



ments. 



Kolh et 



1 2002() used a parallel ring geometry that allowed them to measure 
the normal stress response to shear reversal in concentrated suspensions, in addition to 
measuring the shear stress behavior, and found a transient response in both the normal 
and the shear stresses when the shear was restarted in the opposite direction. Moreover, 
the absolute value of both the normal and the shear stresses changed at the very in- 
stant of flow reversal, which means that the fore-aft asymmetry in the microstructure 
alone is not enough to explain the observed response in the stress upon shear reversal, 
but that non-hydrodynamic forces must also have been acting on the system, either in 
the form of repulsion forces or of rough contacts between spheres. Even more compli- 
cated shear stress responses, including shear-induced ordering, has been recently reported 
at large concentrations {<f> > 50%), probably corresponding to a regime in which non- 



hydrodynamic interactions dominate the behavior of the system l|Voltz et 



\mM. The 
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first direct observations of tlie microscopic structur e developed by dilute suspension s 



Husband fc Gadala-Maria 



(19871), 



{4> = 1% — 5%) undergoing shear were presented by 
who measured in a Couette device the relative distribution of spheres centers in the plane 
of shear, and then by averaging over many realizations, found an anisotropic but fore-aft 
symmetric distribution of close particles. The anisotropy was attributed to the presence 
of pairs of spheres rotating aro und each other forming perma nent doublets. On the other 



hand, in similar experiments. 



Parsi fc Gadala-Marial l|l987ll showed that concentrated 



suspensions (0 — 40% — 50%) do exhibit fore-aft asymmetry, with a larger probability 
of finding pairs of spheres oriented on the approaching side of the reference particle, and 
attributed this asymmetric distribution to either the intrinsic roughness of the spheres or 
to the presence of a n on-hydrodynamic repulsive force between particles. More recently. 



R.amnall et al. 



[ 1997|) used a substantially improved flow-visualization technique to mea- 
sure the pair distribution fimction of dilute suspensions (0 = 5% — 15%) undergoing 
simple shear flow in a shear tank apparatus, and showed that, contrary to the results of 



Husband fc Gadala-Maria 



1 19871) . there is a depletion of permanent doublets moving in 
the region of closed streamlines, and that even for concentrations as small as 5% the dis 
tribution is fore-aft asymmetric. Using the surface roughness model of 



da Gunha fc 



Rampall et al. 



Iin£ 



( 1997 ) 



( 1996|) . and assuming that no particles formed permanent doublets 
were able to reproduce the qualitative trends in the pair distribution function, but the 
predicted depletion of spheres in the regions aligned with the flow was much larger than 
that observed. 

In view of the contradictory results outlined before, it is clear that numerical simu- 
lations, specifically Stokesian dynamics wh ich are well suited for studying low- Reynolds 



number flows of suspensions ((Brady 



200 1|) . offer an important complement to experi- 



ments, in that they can provide detailed, microscopic information that is not accessible 
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via currently available experi mental techniqu e s. In their original work describing their 



Stokesian dynamics method, 



Bossis fc Bradvl (|l98^ showed that the pair distribution 



function of unbounded suspensions undergoing simple shear flow had an angular depen- 
dence, with the microstructure being no longer fore-aft symmetric, and that very few 
particles were oriented in the receding side of the reference sphere. In paper I, we also 
showed such a break in the fore-aft symmetry in the presence of large non-hydrodynamic 
forces acting between the spheres but, for sufficiently small repulsion forces, we found 
that, although anisotropic, the pair distribution function becomes fore-aft symmetric, as 
expected for purely hydrodynamic interactions. To our knowledge, however, as yet no 
systematic numerical investigation has been made of the velocity fluctuations in sheared 
suspensions and their dependence on the underlying microstructure of the suspension. 

It is the purpose of this paper to investigate the velocity fluctuations present in a macro- 
scopically homogeneous, unbounded suspension of neutrally buoyant, non-Brownian sphe- 
res subject to a simple shear flow in the limit of vanishingly small Reynolds numbers using 
Stokesian dynamics, and their dependence on the microstructure developed by the sus- 
pensions. First, we shall focus on the anisotropic, but fore-aft symmetric, distribution of 
close pairs observed in dilute suspensions, and show how it can be accurate ly described as- 



sumin g the absence of permanent doublets of spheres, as first suggested by 



Rampall et al 



1 1997J). We shall also point out that, although the use of a non-hydrodynamic interpar- 
ticle force of extremely short range will yield symmetric distributions, as was reported 
in paper I, the suspensions develop a fore-aft asymmetry for large enough volume frac- 
tions if the range of that force is kept unchanged. Then, we shall show that the pair 



distribution function (?bg(^) obtained by 



Batchelor fc GreerJ l|l972&j) in the dilute limit, 



accurately describes the microstructure in sheared suspensions, in particular the diver- 
gence of the probability density of finding pairs of spheres nearly touching one another. 
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even though it does not account for the observed depletion of closed pairs. Then, by 

making use of the pair distribution function .930(^)7 we shall compute all the temper- 

at ure components in the dilu te limit by numerically integrating the expressions given 



bv lBatchelor fc GreerJ (|1972oD for the particle velocities of two freely suspended spheres 
interacting only through hydrodynamic forces in the presence of a simple shear flow, 
and then compare the results with those obtained from the numerical simulations. Some 
general properties of the temperature tensor valid for isotropic pair distribution func- 
tions will also be discussed. The velocity fluctuations at larger concentrations show the 
effect of the anisotropic structure developed by the flow in that some symmetries of the 
temperature tensor are lost. We also simulate the velocity fluctuations that would result 
from a random hard-sphere spatial distribution of particles in a simple shear flow, and 
thereby are able to further investigate the effects of the microstructure, both its angular 
and radial dependence, on the temperature tensor. In addition, the numerical simula- 
tions provide a full picture of the velocity fluctuations and to this end we shall present 
the probability density functions for all the linear and angular velocity components at 
three different concentrations, showing a transition from a Gaussian to an Exponential 
and finally to a Stretched Exponential form as the volume fraction is decreased. Finally, 
we shall propose a simple correction to the data reduction analysis of the velocity mea- 
surements in LDV experiments, that only depends on the mean angular rotation of the 
spheres in the vorticity direction, and which substantially improves the interpretation of 
the LDV measurements at low volume fractions. 



2. Simulation method: Stokesian dynamics 

We investigate the behavior of suspensions of non-Brownian particles subject to simple 
shear using the method of Stokesian dynamics. A detailed description of the method is 



given in a review by 



MiviustruLt'w 



velocity fluctuations 



Bradv fc Bossia l|1988f l. and the specifics of our simulations were 



already discussed in paper I, hence only a brief discussion is presented here. The method 
accounts for the hydrodynamic forces between solid spheres undergoing simple shear, 
characterized by a shear rate 7, in the limit of zero Reynolds number. In order to simulate 
the behavior of infinite suspensions, periodic boundary conditions in all directions are 
imposed. The simulated cubic cell contains a fixed number of spheres N, related to the 
volume fraction (p hy (j) — {Aira'^ /3)N/V, where V is the volume of the cell. Interactions 
between particles more than a cell apart are included using the Ewald method. A typical 
simulation consisted of = 64 particles sheared over a period of time t ^ 1007^^, and 
all measurements to be reported in this work are for strains ■^t in excess of 50, when 
the system has reached its steady or fully developed state. The motion of the particles 
was integrated using a constant time step At = 10^'^7^^. The results are averaged over 
Nc '~ 100 different initial configurations, with each initial configuration corresponding 
to a random distribution of non-overlapp ing spheres in the simulatio n cell, using the 
random-phase average method proposed by 



Marchioro fc Acrivoi 



1 20011) . In what follows. 



we shall express all the variables in dimensionless units, using the radius of the spheres 
a as the characteristic length and 7^^ as the characteristic time. 

In a suspension of monodisperse spheres undergoing simple shear the separation be- 
tween spheres may become exceedingly small during two-particle collisions (less than 10^^ 
of their radius), and the effects of surface roughness or small Brownian displacements 
cannot be neglected. Usually, a short-ranged, repulsive force is introduced between the 
spheres to qualitatively model the effect of these non-hydrodynamic interactions, with 
the numerical advantage of preventing any overlaps during close encounters between 
particles. As in paper I, we used the following standard expression for the repulsive 
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interparticle force, 



af3 — — -; zirnr^ap 1 (2.1) 

1 — e '^1^" 



where 67r/ia^7FQ^, with /i being the viscosity of the suspending hquid, is the force exerted 
on sphere a by sphere /3, _fo is a dimensionless coefficient reflecting the magnitude of 
this force, is the characteristic range of the force, e is the distance of closest approach 
between the surfaces of the two spheres divided by a, and e^^ is the unit vector connecting 
their centers pointing from (3 to a. 

The effect of the characteristic range of the interparticle force Tc on the microscopic 
structure of the suspension was discussed in paper I. First, we showed that the mini- 
mum separation reached by colliding spheres, and therefore the flrst peak in the pair 
distribution function, is strongly affected by the range of the interparticle force in that, 
as Tc increases, the minimum separation between neighboring particles also increases. 
Then, we showed that, in general, the presence of a repulsive force breaks the fore-aft 
symmetry of the particle trajectories in a simple shear flow. However, we also showed 
that the symmetry is recovered, for small enough values of the force range, Tc ~ 10""*, 
at least in the sense that no asymmetry was observed in the angular dependence of the 
numerically computed pair distribution function. In this work, we use this small range 
for the interparticle force, Tc ~ 10^''. 

3. Microscopic structure induced by the shear 

The investigation of the microscopic st ructure developed by suspe nsions undergoing 



shear flow followed the pioneering work by 



Batchelor fc GreenI l|19726j) , where an expres- 



sion was derived for the pair distribution function g{r) in the dilute limit. This function 
is related to P(r|ro), the probability of finding a sphere with its center at r given that 
there is a sphere with its center at ro = 0, by P(r|ro) = (/)g(r)(3/47r). Recall that, even 
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a random hard-sphere distribution leads to correlations in the position of any two parti- 
cles due to excluded volume effects, and typically displays a liquid-like microstructure at 
high volume fractions. But, in addition to these excluded volume eff ects, hydrodynamic 



interac tions between spheres lead to surprising results. Specifically, 



Batchelor fc Green 



1 1972fcfl showed that the pair distribution function is an isotropic function of the dis- 
tance between the two spheres, i. e. g{r) — g{r), and that it diverges as r — > 2, which 
means that pairs of particles are substantially more likely to be found near contact in 
a sheared suspension than in a random hard-sphere distribution. In turn, this implies 
a high correlation in the position of the spheres that is not present in a random ha rd- 



sphere configuration. The expression for g{r) derived by 
the dilute approximation is. 



Batchelor fc GreenI l)l9726l) in 



.9bgW 



1 



1 - A{r) 



exp 



3 B{0 - A{0 



(3.1) 



where the mobility functions A and B are functions only o f r. H ere, we shall use the 



expressions for these functions given by 



da 



terval r > 2 into three different regions (see 



^jjnh^^^Hingh J 199^ ), who divided the in- 



da Cunha fc Hindi 



19961) for details on how 



they obtained the expressions for A and B in each region). Specifically: a) within the 

lubrication region 2 < r < 2.01, 

_ (16.3096- 7.1548r) 
A — , 

r 

^ 2 (0.4056L2 + 1.49681L - 1.9108) 
^ r (L2 + 6.04250L -h 6.32549) ' 

where L — — ln(r — 2), 

b) within the intermediate region 2.01 < r < 2.5, 



A 
B 



-4.3833 + 17.7176^ + 14.8204^ - 92.447l4 - 46.3151-"'^ 



-3.1918 + 12.364li + 11.4615^ 



65.2926- 



36.4909- 



232.2304- 



154.8074- 
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and c) within the far- field region r > 2.5, 



25 


35 


125 


102 


12.5 










^11 
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25 


36 \ 




rp8 











3 \ r' 

To be precise, these resuhs for the pair distribution function gsci'''): f^pply only to 
particles that are initially far from each other, and therefore spatially uncorrelated. But, 
in the limit of purely hydrodynamic interactions and very dilute suspensions (no three- 
particle interactions), there is a region of close trajectories in r-sp ace, where pairs of 



particl es remain correlated at all times forming permanent doublets l|Batchelor fc Green 



1972qI) . In this case, accounting for the probability distribution of particles forming per- 
manent doublets would require the knowledge of the initial distribution of the particles. 
On the other hand, the presence of any non-hydrodynamic interaction, such as roughness 
or repulsion forces, or the existence of three-particle collisions, would generate a transfer 
of particles across the streamlines and therefore remove the need for specifying the initial 
distribution of spheres. Even so, to obtain the probability distribution would still require 
the full knowledge of the transfer process (the combination of three-particle collisions and 
non-hydrodynamic interactions), and the solution of the corresponding boundary value 
problem. 

This implies that the pair distribution function may actually be anisotropic simply due 
to the distribution of pairs forming permanent doublets in that, although the distribution 
of particles outside the regio n of closed streamlin e s does not have an angular structure in 



the dilute limit, as shown by 



Batchelor fc GreenI l|l972a) , the distribution of particles in 



the region of closed streamlines, which extends to r oo, may actually render the com- 
plete pair distribution function anisotropic. In fact, in paper I we showed that, although 
for exceedingly short ranged repulsion forces, ~ 10^"*, the pair distribution function 
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180 135 90 45 



approaching q ^^j^g receding 

Figure 1. Normalized angular distribution function gR{6) for pairs of particles. The distance 
between the pairs lies in the range 2 < r < 2.01 (R = 2.01). Different curves are for different 
volume fractions. All simulations were performed with A'^ — 64, Nc = 100, Fo — 1 and rc = 10~*. 
The solid line is obtained using the pair distribution function given in Eg 13. II for the region out- 
side the closed streamlines and assuming zero probability of finding a pair forming a permanent 
doublet (see text for a more detailed explanation). 

for close particles recovered its expected fore-aft symmetry, it remained anisotropic. In 
figure n we present gR{0), the angular dependence of the pair distribution function of 
pairs closer than a certain distance R, as defined in paper I, for different concentrations of 
the suspension (as mentioned in section |21 numerical results are for a range Tf, = 10~* 
of the interparticle force, which is the smallest value of the force range simulated in paper 
I). It can be seen that, even for this exceedingly short ranged interparticle force, fore-aft 
symmetry is broken at large enough concentrations, and that the distributions show a 
larger number of pairs oriented on the approaching side of the reference sphere than on 
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x2 




^ Closed Streamlines Boundary 



xl 



Figure 2. Schematic representation of the computa tion of Qnid) using the regi ons of open and 



1972a). The sohd hue 



closed streamhnes for a pair of interacting spheres (iBatchelor fc GreenI 
represents the intersection of the plane of shear xi — X2 with the surface bounding the region of 
closed trajectories, which can be formed by rotating this curve about the X2 axis together with 
its mirror image obtained by reflection about the xi axis. The gray circle of radius a represents 
the reference sphere and the circle of radius 2a encloses the excluded volume. Then, a particle 
with its center located inside the sphere of radius R forms, with the reference sphere, a pair 
that is closer than R, and is included in gR{9). For a given angle 9, the distribution of pairs 
closer than R is calculated by integrating the pair distribution function given in Eq. 13.11 only 
in the shaded region, which corresponds to open trajectories only, because the region of closed 
trajectories is considered to have a negligible effect due to the depletion of permanent doublets. 

the receding one. More importantly, the angular distribution function seems to reach, in 
the dilute hmit, an asymptotic distribution which is anisotropic and shows a depletion of 
pairs oriented close to the flow direction. This suggests a depletion of permanent doublets, 
which spend more time nearly aligned along the direction of the flow, and seems to indi- 
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Figure 3. Angular distribution function gR{6) for pairs of particles corresponding to i? = 2.01 
and R = 2.1. The volume fraction is (f) = 5%. The solid line is obtained using the pair distribution 
function gaGir) for the region of open trajectories and assuming zero probability inside the region 
of closed trajectories, (see text and figure |21 for a more detailed explanation). 



cate that, as speculated by 



Rampall et al. 



(I993), any mechanism forcing particles into 



the region of closed streamlines is small compared to the effect of the non-hydrodynamic 
forces which eliminates particles from this region, and ultimately leaves only a negligible 
number of pairs forming permanent doublets. In this case, the pair distribution function 
in the dilute limit should be the combination of guci'i') for the region of open trajectories 
and a zero probability inside the region of closed streamlines. We show schematically, in 
figure 13 how we can then approximate the angular dependence of the pair distribution 
function of pairs closer than a certain distance i?, using the expression for the surface 
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separating the regions of open and closed trajectories, ( 



id A. Avi ivub^ 



Batchelor fc Green 



Clearly, the surface is axisymmetric with X2 as the symmetry axis. (Here, and in what 
follows, the Cartesian axis 1 lies along the direction of the mean flow, 2 is perpendicular 
to 1 along the plane of shear, and 3 is the vorticity axis.) 

In figure ^ we compare this approximation to the angular dependence of the pair 
distribution function of close pairs (i? = 2.01) with the numerical results, and find a very 
good agreement for concentrations smaller than 5%. Moreover, in figure^ we show that 
this approximation accurately describes the anisotropy found for the pair distribution 
function of pairs with an order of magnitude larger range of separations, i.e. R = 2.1, thus 
validating the assumption of a complete depletion of pairs forming permanent doublets. 

Finally, in figure 01 we present the radial dependence of the pair distribution function, 
i.e. the pair distribution function integrated over both spherical angles, as obtained in the 
numerical simulations at different particle concentrations. We also compare the numerical 
results with the pair distribution function given in Eg. 13. II and find that, although .gBol'") 
does not account for the effect of the closed trajectories, in particular the observed 
depletion of permanent doublets, it both follows the simulation results fairly accurately 
over a wide range of r, as well as captures the substantial increase in the probability of 
finding pairs of particles near contact. On the other hand, 5bg('') seems to overestimate 
the asymptotic distribution in the dilute limit, and is unable to take into account the 
effects of the depletion of permanent doublets, which would ultimately lead to g{r) — 
for r smaller than the minimum possible separation between approaching spheres in the 
region of open trajectories, rmin ~ 4 x 10^^. But, for distances between the spheres that 
are not too small, r > 2.001, Eg. 13.11 captures the divergent trend of g{r) as r — > 2, 
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Figure 4. Pair distribution function integrated over all possible angular orientations. Differ- 
ent symbols correspond to the results obtained in numerical simulations for different volume 
fractions, (ft = 0.25, 0.15, 0.10, 0.05 and 0.03. The pair distribution functions are constructed 
from a histogram of the distance between all pairs of particles, averaged over time and over 
different realizations (the smallest size of the bins in the histogram is Ar = 0.005). The solid 
line corresponds to the pair distribution function given in Eq. 13.11 

which as we shall see, allows us to obtain a reasonably accurate estimate of the velocity 
fluctuations in the dilute limit. 



4. Velocity fluctuations 



Following Batchelor's notation IIBatchelo: 



1^ 



matrix of the velocity fluctuations (^van Kampen 



the t emperature tensor (the covariance 



1987j) ) can be written as, 



T 



dCN 5v,{r)5vj{r)P[CN\r), 



(4.1) 
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where Svi{r) is the fluctuation in the velocity component Vi for a particle located at r 

when the configuration of the surrounding spheres is given by Cn, with P(Cjv|t') being the 
probability of such an event. Prom its definition, it is clear that the temperature tensor is 
symmetric = Tjj. In addition, in simple shear flows there exists an inversion symmetry 
in the vorticity direction in that, a given configuration and its counterpart in which X3 is 
changed by —X3 are equally probable, and therefore we have that T13 = T23 = 0, for any 
volume fraction. We can simplify the temperature tensor even further by decomposing 
the simple shear flow into a solid body rotating flow, which does not contribute to the 
velocity fluctuations irrespective of the concentration of particles, and a purely straining 
motion. The latter is symmetric in xi and X2 and therefore, for any particular velocity 
fluctuation, say in the 1 direction, in a conflguration C_sr of particles surrounding the 
reference sphere, the same fluctuation but in the 2 direction would be obtained by a 
configuration Cj^ in which all the particle positions in Cjv are transformed according 
to Xi X2- Then, it clearly follows that Tn = T22, depending only on whether the 
configurational probability density P{CN\r) has the same symmetry, i.e. it is invariant 
under the transformation xi ^ X2- Moreover, since a configuration C'Jj- in which the 
particle positions in Cj^ are transformed according to xi —xi would give the negative 
of the previous velocity fluctuation, and similarly for fluctuations in the 2 direction, it 
is clear that T12 = 0. Thus, the temperature tensor should be diagonal, as long as the 
probability density of particle configurations is invariant under those changes, i.e. as long 
as P{CN\r) remains invariant under the transformations Xi <-> —Xi and X2 <-»■ —X2- 

In summary, for any concentration and a symmetry preserving configurational proba- 
bility P{CN\r), we have that the off-diagonals terms of the temperature tensor are null 
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and that the temperatures in the plane of shear are equal, 

=0 i^ j 
Til = T22 
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(4.2) 
(4.3) 



In the dilute limit, the fluctuations in the velocity come from two-particle interactions, 
and from the far-field form of these interactions it can be shown that any component 
of the temperature tensor, of the form 6vi6vj, decays faster than \/r^ and therefore, 
its average value can be directly computed by averaging the hydrodynamic interaction 
between a pair of spheres over all possible configurations, 



^ J dr Sv,Svjg{r) 



= 4) Uj, 



(4.4) 



which gives a linear dependence of the temperature components on the volume fraction. 



For two freely-moving spheres in a simple shear flow, the velocity fluctuation of a 
sphere induced by a sec ond sphere the center of which is located at r is given by 



I da Cunha fc Hinch 



199a): 



Svi — xi — X2 — exi — 2^^^ 
6v2 = ex2 - ^Bxi 
Svs = ex3, 



(4.5) 
(4.6) 
(4.7) 



where e = XiX2{B — A)/r'^. Using these equations, it can be easily shown that for a 
pair distribution function which depends only on r, the temperature tensor is not only 
diagonal, but that T33, its component in the vorticity direction, is smaller than the 
temperature in the plane of shear. 



Til — T22 > T- 



33- 



(4.8) 



The exact temperature values will depend in general on the pair distribution function. 
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til ~ t22 


^33 






^33 


Random Hard Sphere (gusir)) 


0.3f57 


O.OSff 


Simple Shear Flow (5bg(»')) 


0.4637 


0.1031 


lubrication 


0.0040 


0.0006 


lubrication 


0.0896 


0.0117 


intermediate 


0.0930 


0.0f75 


intermediate 


0.i53f 


0.0279 


far-field 


0.2f87 


0.0630 


far-field 


0.2210 


0.0635 



Table 1. Temperature tensor in the dilute limit, computed using Egs. 14.414.71 for two dif- 
ferent pair distribu t ion fun ctions, one for a hard-sphere distribution, and the other given by 



Batchelor &l Green 



( 1972a) for a simple shear flow in the dilute limit, Eq. 13.11 Also given is 
the contribution to the velocity fluctuations coming from the three different regions in which 
the mobility functions A and B are divided, i.e. the lubrication region (2 < r < 2.01), the 
intermediate region (2.01 < r < 2.5) and the far-field region (2.5 < r). 



In Tabled we present the diagonal terms of the temperature tensor in the dilute limit, 
obtained from the numerical integration of Eq. 14.41 for two different isotropic pair dis- 
tribution fimctions, corresponding to a random distribution of hard spheres, gnsi''') = Ij 
and to Batchelor fc Green's result f or a suspension in a simple shear flow, ^bg (j") given by 



Ea. l3.1l f Batchelor fc Greenlll972M) . Note that, although in our numerical simulations we 
observe a depletion of permanent doublets in the dilute limit, we first neglect this effect, 
and compute the temperature terms by numerically calculating the integrals in Eg 14.41 
using gBoC^), for all possible angular orientations. As we shall see, for dilute suspensions, 
this provides a satisfactory approximation to the temperature tensor computed from 
our numerical simulations. The contribution to the integral of each region in which the 
mobility functions A and B are divided, i.e. the lubrication, intermediate and far-field 
regions, is also given in Table ^ As expected, the far-field contribution is practically 
identical in both cases, since in fact 5bg(^) asymptotically approaches gnsif) for r — > oo. 
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Figure 5. Diagonal components of the temperature tensor as a function of the volume frac- 
tion, obtained from the numerical Stokesian dynamics simulations. The solid and dashed lines 
correspond to the dilute limit calculation for Tii-22 = iii;22 and Tss — (j) ^33, respectively. 
The computed values of tii;22 and 133 are given in the second part of tabled The discrepancy 
between theory and numerical results is about 25%. 

On the other hand, the contribution of the lubrication region to the velocity fluctuations 
is at least an order of magnitude larger if computed using (7bg('') because, in that case, 
the probability of finding two nearly touching spheres is substantially larger than in the 
hard-sphere case. On the other hand, the anisotropy ratio in the dilute limit is similar in 
both cases, Tii/r33 ^ 4. 

In figure |S1 we present the diagonal terms of the temperature tensor as a function of 
the volume fraction, obtained in a simple shear flow by means of Stokesian dynamics 
simulations. The temperature components Tn and T22 converge to a common curve in 
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Time 



Figure 6. Velocity autocorrelation function in the dilute limit in both transverse directions. 
Symbols correspond to the results obtained using the SD method for a volume fraction — 0.03. 
Solid lines correspond to the numerical computation of the velocity autocorrelation function 
using two-particle interactions only. 

the dilute limit, whieh is eonsistent with the existence of an isotropic pair distribution 
function and indicates that the effect of the particle-depleted region of closed streamlines 
is not measurable. In addition, the decay of the velocity fluctuations follows the dilute 
limit scaling given by Eq. 14.41 viz. that Tij is proportional to </>, even for surprisingly 
high volume fractions. On the other hand, at larger concentrations we see that the Tn 
and T22 curves separate from each other, which is evidence of the structure developed by 
the suspension at high concentrations. In fact, in figure ^ we showed that, although we 
have a very short-ranged interparticle force, Vc — 10^''', the pair distribution function has 
fore-aft symmetry in the dilute limit, for larger concentrations this symmetry no longer 
holds. 
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The agreement observed in figure El between the suspension temperature computed 

in the dilute limit using .(7bg(^) firid via the numerical simulations, shows that, as ex- 
pected, the dominant contribution to the velocity fluctuations arises from two-particle 
hydrodynamic interactions. Furthermore, in paper I, we argued that in the dilute limit, 
the whole velocity autocorrelation function converges to an asymptotic function domi- 
nated by two-particle encounters. To further investigate the dilute limit behavior of the 
velocity fluctuations, we therefore computed the velocity autocorrelation function on the 
basis only o two-particle hydrodynamic interactions, using Eqs. 14.514.71 In order to do 
that, we simulated a large number of two-particle encounters {N^ — 2 x 10^) between a 
test sphere, initially located at the origin, and an incoming particle, initially far away 
from the test sphere. The exact position of the incoming particle was chosen randomly 
in the region {~/S.xi < < —/S.xi 4- i?, < < i?, < a;^ < i?), where R is the size 
of the cross-section for the incoming particles considered in the calculations. Note that, 
in order for an encounter to induce a significant velocity fluctuation in the test sphere, 
both spheres should come rea sonably close to ea ch other at some point during their mo- 



tion, and we thus use R = 5 l|Wang et 



1996l) (Axi was set to 20R). The probability 



distribution used to generate the initial conditions was uniform in xi and X3, and in 
the shear direction we used a probability distribution proportional to the incoming flux 
of particles in simple shear flow, that is p{x2) oc X2. The cross-section of the region of 
closed streamlines, perpendicular to the flow direction, is so small at Axi — 20R, that 
we did not observe any closed trajectory after simulating Nc encounters. The motion of 
both spheres was then computed, using Eas l4.5lC7l and a time step At = 10^^, until 
the incoming particle reached the point that was symmetric with respect to its initial 
position, that is {—xi, Xj, Xg). Finally, the velocity autocorrelation function was com- 
puted by averaging over all the simulated trajectories, after splitting each one of them 
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into intervals of time T = 10. The results are shown in figure El where the numerically 

computed velocity autocorrelation functions in both transverse directions is compared 
to the results obtained by means of Stokesian dynamics simulations (already presented 
in paper I). An excellent agreement is obtained, which confirms that the velocity auto- 
correlation functions in both transverse directions reach corresponding asymptotic forms 
in the dilute limit. It also demonstrates that, as was first suggested in paper I, the fact 
that both autocorrelation functions become negative at times t ~ 1 is due to the anti- 
correlated motion performed by the spheres during binary collisions, i.e. the transverse 
velocities of the spheres involved in a binary collision are reversed at the instant at which 
the incoming particle goes from the approaching to the receding side of the reference 
sphere. 

In order to investigate the effects of the steady state microstructure developed by the 
sheared suspensions, we performed a second type of numerical computations, in which the 
velocity fluctuations were calculated for a random hard-sphere distribution of particles 
subject to simple shear flow. That is, we first generated a random distribution of spheres 
at the desired concentration, and then, using the Stokesian dynamics code, we calculated 
the instantaneous velocity of all the spheres in the presence of a simple shear flow. We 
then averaged the results over many different realizations of the random hard-sphere 
distribution. Typically, the number of particles in each realization was the same as in the 
dynamic simulations, but the number of configurations was increased to Nc — 256. We 
refer to the randomly generated hard-sphere particle distribution as Hard Sphere (HS) 
distribution in contrast to the Shear Flow (SF) distribution which refers to the particle 
distribution which is attained asymptotically after the suspension has been sheared in a 
simple shear flow for strains in excess of 50. 

In figure we plot the diagonal components of the temperature tensor, obtained for 
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Figure 7. Diagonal components of tlie temperature tensor as a function of the volume fraction, 
as obtained for a random distribution of liard splieres. We assume tliat tiie splieres are subject to 
a simple shear flow and compute their velocity using the Stokesian dynamics code. The dashed 
and dotted lines correspond to the dilute limit calculation Tij — (j> tij and tij given in the first 
part of table Q 

HS distributions in simple shear flow. We can see that, as a result of the isotropic spatial 
distribution of particles, Tn equals T22 for all values of the volume fraction (within 3%). 
Also, r33, the temperature in the vorticity direction is smaller than that in the plane 
of shear, as predicted in the dilute limit. The agreement between theory and numerical 
results is excellent in this case, with a discrepancy less than 10% for the lowest volume 
fraction. 

As mentioned earlier, the off-diagonal terms of the temperature tensor are expected to 
be zero for an isotropic pair distribution function. Moreover, if the only broken symmetry 
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Figure 8. Off-diagonal terms of the temperature tensor as a function of the volume fraction 
is the fore-aft symmetry, the only term that may dilfcr from zero is T12, which should 
actually be negative if particles are depleted in the receding side of the reference sphere, 
as is observed in the numerical and experimental work. The numerical results agree 
completely with this analysis, as shown in figure |H1 in that, for all volume fractions, all 
the off-diagonal components vanish for the HS distributions, as well as Tia and T23 for 
the SF distributions. Moreover, the only correlation present in the velocity fluctuations 
is given by T12, which becomes different from zero only as the concentration increases 
and a fore-aft asymmetry is developed by the SF distributions. 

A completely analogous analysis to the one presented above for the linear veloc- 
ity fluctuations gives very similar results for the fluctuations in the angular velocity, 
Qij ~ (SuJiScuj). The off-diagonal terms are zero (flij — 0, i ^ j) and Qn ~ Q22, for 
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Figure 9. Diagonal components of the angular velocity autocorrelation tensor as a function of 
the volume fraction, obtained from the Stokesian dynamics simulations. The solid and dashed 
lines correspond to the dilute limit calculation for Q,\i-22 ~ nii;22 and ^33 — (f> Cl-n, respec- 
tively. Here the black and gray lines correspond to the dilute limit calculations using SF and 
HS distributions respectively, and the corresponding values of Cli\-22 and Has are given in table 
|5| The discrepancy between theory and numerical results is close to a factor 2. The results are 
bounded between the calculations using SF and HS distributions. 

all concentrations, as long as the distribution of spheres has the symmetries discussed 
previously when we analyzed the properties of the temperature tensor. 
In the dilute limit, the fluctuations can be written as, 

3 



f^ij — dr SuJiSujj (3(/)/47r) 5(r) = (p 



— j dr 5uJi5ujjg{r) 



= 4) 



(4.9) 
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fill = fi22 = fisa/S 




fill = fi22 = fiss/S 


Random Hard Sphere {gus{r)) 


0.01064 


Simple Shear Flow (s'bg('')) 


0.02260 


lubrication 


0.00033 


lubrication 


0.00869 


intermediate 


0.004438 


intermediate 


0.007915 


far-field 


0.00587 


far-field 


0.00600 



Table 2. Angular velocity fiuctuations in the dilute limit, computed for two different assumed 
pair distribution functions, one corresponding to a random distribution of hard spheres, and 
the other using i7bg('"). The contribution to the velocity fluctuations coming from the three 
different regions in 5bg('") and C{r) - c.f Eg. 14.131 - are noted separately, i.e. the lubrication 
region (2 < r < 2.01), the intermediate region (2.01 < r < 2.5) and the far-field region (2.5 < r). 



SuJ2 



2" -2 
-C- 



1 — X2 



and for two freely-moving spheres we have that (|Batchelor fc GreerJll972(7ll : 

OLOl = --C 5—, 



(4.10) 
(4.11) 
(4.12) 



where C is a dimensionless function of r only. Using Eqs. 14. 914.121 and assuming an 
isotropic pair distribution function, it can be shown that, fJaa — 5 fJn = 5 ^22 for 
any radial dependence of the pair distribution fun ction. Finally, u s ing th e far-field and 



Kim fc Karrilal ((l99l|), and a linear 



the lubrication approximations of C(r) given by 
interpolation in the intermediate region, 

^ly^l ~ 6. 32549+6. 0425L+L^ ? < /.Ui 

C^{r)^^^,Cl{r)^^^^Cf{r) for 2.0K r < 2.5 
Cf{r) = J^-^+^^ + ^ + ^ for 2.5<r 



C{r) = I 



(4.13) 



we obtain the results presented in table [5] for ^li 
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<l> 

Figure 10. Diagonal components of the angular velocity fluctuations as a function of the volume 
fraction, obtained for HS distributions, subject to a simple shear flow, and computed using the 
Stokesian dynamics code. The solid and dashed lines correspond to the dilute hmit calculation 
— (j) ^ij with being given in the first part of table|5| The difference between theory and 
numerical simulations is less than 20% at the lowest concentrations. 

In figures |51 and ^1 we present the velocity fluctuations in the angular velocity ob- 
tained for the SF and the HS distributions of particles. These fluctuations seem to be 
very sensitive to the microstructure developed by the suspension because, even at very 
low concentrations, a difference between fJn and f222 can be observed in the SF case, 
contrary to what happens in the HS case where they coincide for all concentrations, as 
expected from our previous discussion. Also, the discrepancy between the theoretical and 
the numerical values is large for the SF distribution (a factor ^ 2), while, in the HS case, 
there is good agreement with the theory. Not even the linear behavior on (p in the dilute 
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Figure 11. OfF-diagonal terms of the angular velocity autocorrelation tensor as a function of 

the volume fraction 

limit seems to have been reached in the SF case at low concentrations, in contrast to 
what is observed for the HS distributions of spheres. However, we also show in figure El 
that, the calculations using gHs(^) provide a lower bound to the angular velocity fluc- 
tuations. Qualitatively, this behavior can be understood from the observation that the 
main difference between the results obtained for HS and SF distributions comes from 
the contribution of the lubrication region to the fluctuations, which is negligible in the 
HS case. Thus, a lower bound to the velocity fluctuations in the absence of permanent 
doublets can be estimated roughly using the HS distributions, which have a negligible 
contribution from the lubrication region. 

Finally, we present the off-diagonal terms of the angular velocity autocorrelation ten- 
sor. As predicted, all terms are negligible in the dilute limit, and only the microstructure 
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Figure 12. Probability density functions for the three components of the linear velocity fluctu- 
ations and the three components of the angular velocity fluctuations, at very low concentrations 
(f) ~ 0.01. All p.d.f.'s for the indicated variables, i;2, ^s, tiJi, ^1^2 and Slus, are rescaled using the 
corresponding temperatures so that the standard deviation is unity (for example ^ = Svi / \/Tii, 
and similarly for the other variables). 

developed by the SF distributions leads to a correlation ^ 0, which in fact is in agree- 
ment vi^ith the result, referred to earlier, that ti2 7^ 0, and with the observed depletion 
of pairs of particles oriented on the receding side of the interaction (see figure QJ. 

In addition to the temperature values, the numerical simulations provide us with 
greater detail about the velocity fluctuations. In fact, we can obtain the full proba- 
bility distribution function for the fluctuations in both the linear and angular velocity 
components, calculated from a histogram of the particle velocities, averaged both over 
Nc different realizations and in time. 
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Figure 13. Probability density functions for the three components of the linear velocity fluc- 
tuations and the two components of the angular velocity fluctuations in the plane of shear. 
The vorticity component of the angular velocity fluctuations is not symmetric as it can be ob- 
served in flgure \Wl The volume fraction is = 0.10. AU p.d.f.'s for the indicated variables, 
5vi,V2,V3,u;i,u>2 and Suj^, are rescaled using the corresponding temperatures so that the stan- 
dard deviation is unity (for example ^ — Svi/^/Tu, and similarly for the other variables). 
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UJi 
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0.45 


Vi 


2.54 


UJl 


16.1 


0.10 


Vl 


14.1 


UJl 


39.2 


0.01 


Vl 


18.1 


UJl 


22.8 


P = 2 


V2 


1.85 


LJ2 


17.8 


P = l 


V2 


10.9 


UJ2 


36.0 


13 = 1/4: 


V2 


18.1 


UJ2 


21.7 




V3 


4.54 








Vg 


21.1 








V3 


24.8 







Table 3. Fitted values for the probability density function of all the linear and angular 

velocity components. 
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Figure 14. Probability density functions for the three components of the hnear velocity fluc- 
tuations and the two components of the angular velocity fluctuations in the plane of shear. 
The vorticity component of the angular velocity fluctuations is not symmetric as it can be ob- 
served in flgure [T^ The volume fraction is = 0.45. All p.d.f.'s for the indicated variables, 
5v\,V2,V2.,uJ\,uj2 and 5lj3, are rescaled using the corresponding temperatures so that the stan- 
dard deviation is unity (for example ^ = 5wi/\/Tn, and similarly for the other variables). 

In figurcs [T^ll3l and ll4l wc show the normalized p.d.f of the linear and angular velocity 
fluctuations in all directions, and for three different volume fractions. (Note that, for the 
velocity component in the direction of the shear, we have subtracted the contribution of 
the external velocity field at the center of the particle, that is 5vi = ii — X2.) Different 
functional forms are observed as the concentration decreases. A first transition, from 
Gaussian to Exponential distributions occurs when the concentration is decreased from 
(j) = 0.45 to (j) = 0.10, as already presented in paper I. A second transition, from Expo- 
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Figure 15. Probability density functions of the angular velocity in the vorticity direction, 
P{u!3), for five different volume fractions (ji = 0.45, 0.35, 0.25, 0.15 and 0.05. 

nential to a Stretched Exponential with exponent ~ 1/4 occurs when the concentration 
decreases even further down to (j> = 0.01. All the numerical data were fitted using expo- 
nential distributions of the form P(S,) oc ea;p(— ^"), with a = 2 for large concentrations 
{(p — 0.45), a = 1 for intermediate values of the volume fraction (0 — 0.10), and a = 1/4 
for very low concentrations {(f> = 0.01). In paper I, we discussed the first transition, from 
Gaussian to Exponential distributions, by analogy with turbulent flows, where one ob- 
serves this type of transition in the p.d.f of the velocity differences, the temperature and 
other passive scalars. The second transition presented here is in accordance with that 
analogy, in that, with decreasing concentrations, the exponent a of the distribution also 
decreases, implying intermittency, in th at the probability of rare events is much larger 
than expected from Gaussian statistics 
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Figure 16. Mean value of the angular velocity in the vorticity direction {ui^} as a function of 
the volume fraction (j) of the suspension. Open symbols corresponds to the same mean angular 
velocity but for the HS distributions, calculated using the Stokesian dynamics code. 

Finally, the vorticity component of the angular velocity, ^3, has a mean value dif- 
ferent from zero, due to the shear flow. In the Stokes limit, by decomposing the linear 
shear into a purely rotational flow and a purely stra ining flow, it is easy to show that 



1992), which is therefore the 



the angular velocity of a single sphere is fl^ ~ ~l/2 l|Lea]| 
expected average value in the dilute limit. For larger concentrations, however, hydrody- 
namic interactions between spheres should be taken into account. But, using the same 
superposition of flows, and due to the reflection symmetry of the purely straining flow, it 
can be shown that the average remains constant, and equal to fis, even in the presence 
of other spheres, as long as the distribution of spheres is isotropic. In figure [T31 we show 
the distribution of angular velocities in the vorticity direction, for (j) — 0.45, 0.35, 0.25, 
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0.15 and 0.05, while in figureEl we show the mean angular velocity as a function of the 

volume fraction. As can be seen, (^3) decreases from —1/2 down to —0.58 at (p = 0.35 and 
then increases to —0.55 at </> = 0.45. (However, note that the shift in (wa) with respect to 
= —1/2 is always smaller than the width of the distribution a^^, and therefore that 
the fluctuations are larger than the shift in the average value.) We also show the results 
obtained for the HS distributions, calculating (wa) using Stokesian dynamics, where it 
can be seen that the mean angular velocity remains equal to fis = — 1/2 for all concen- 
trations. We can conclude therefore, that the deviation in the mean angular velocity is 
due to the anisotropy developed by the suspension in simple shear flows. Moreover, as 
shown in figure 3 of paper I, the angular dependence of the pair distribution function 
for close spheres shows a larger probability for pairs oriented at angles 45° < 9 < 135°, 
which is consistent with an increase in the angular speed of the spheres. 

4.1. A note on LDV measurements 

It is known that, due to the spatial but random distribution of the scattering sites within 
the spheres, LDV measurements contain spurious contributions to the linear particle ve- 
locity, resulting from the rotation of the particles, which are invariably neglected. This 
would appear to be permissible for the case of the temperature measurements, given that 
most of these spurious contributions average to zero and do not affect the variance of the 
velocity fluctuations because th e location of the sca ttering sites is uncorrelated from one 
particle to another. In addition. 



Lvon fc Lea 



( 19981) estimated the contribution from the 



mean particle rotation to be one order of magnitude lower than the velocity fluctuations 
resu lting from interpartic le interactions, and based upon this argument neglected its ef- 



fect. 



Shaplev et al. 



(|2002() explicitly computed the contribution of the average rotation of 
the spheres to the measured velocity fluctuations, but also concluded that its magnitude 
was negligible compared to the fluctuations resulting from collisions between particles. 
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Figure 17. Particle velocity fluctuations in the mean flow direction, Tn, as computed directly 
from the numerical simulations (solid circles), as if they were measured using the LDV technique 
(open circles) and after correcting the LDV measurements to account for the average rotation 
in the vorticity direction (open triangles). 

However, the spurious contributions to the measured velocity fluctuations originating 
from the mean angular velocity in the vorticity direction are independent of concentra- 
tion in the dilute limit, and moreover, we have shown earlier that, again in the dilute 
limit, the angular velocity fluctuations are proportional to the volume fraction. Thus, 
it is clear that the spurious contribution to the measured velocity fluctuations due to 
the angular rotation of the spheres eventually becomes important, and even dominant, 
at low enough concentrations. For example, if the scattering sites were distributed uni- 
formly inside the spheres, which rotate with mean angular velocity {0J3), the measured 
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SD of the velocity fluctuations in the direction of the flow, T^^ , can be written as, 



where the second term on the right hand side corresponds to the above mentioned spu- 
rious contribution to the measured velocity fluctuations due to the rotation of the par- 
ticles, and is usually neglected. But since, in the dilute limit, Tn + (1/5) (5122 + ^^33) = 
[til + (l/5)(fi22 + f^33)]0 (c.f. Eq.OIl, while (0)3) = Viz = -1/2, it is clear that, at low 
concentrations, the mean angular rotation of the spheres dominates the contribution to 
Tii^^ , hence, subtracting ^13/6 — 1/20 is a correction needed in this limit. To illustrate 
this result we compare, in figure [T7I the velocity fluctuations in the direction of the flow 
as would have been obtained from the LDV measurements, Tj"^^ , and that from the 
same measurements but corrected by the average rotation, Tf"^^ — with the real 

temperature Tn. Surprisingly, this simple correction to the analysis of the data remains 
important even at concentrations as large as 20%, or even larger, and in fact the cor- 
rected values for T^^^ stay very close to the true temperature Tn over the whole range 
of concentrations investigated. 

5. Summary 

The velocity fluctuations that occur when a simple shear is imposed in a macroscop- 
ically homogeneous suspension of neutrally buoyant, non-Brownian spheres, and their 
dependence on the microstructure developed by the suspensions, were investigated in 
the limit of vanishingly small Reynolds numbers by means of Stokesian dynamics simu- 
lations. These simulations account for the hydrodynamic interactions between spheres, 
and also include a short-range repulsion force that qualitatively models the effects of 
surface roughness and Brownian forces. We simulated the evolution of a large number of 
independent initial hard-sphere random distributions for strains 7i ~ 100, which in our 




(4.14) 



Micro structure and velocity fluctuations 39 
previous work, proved sufRcie ntly long to allow u s to study the system in the asymptotic, 



fully developed steady state IjDrazer et al 



2002) 



We first discussed the angular structure developed by the suspension undergoing sim- 
ple shear, and showed that, even for exceedingly short ranged interparticle forces, the 
distribution of particles is fore-aft asymmetric at large concentrations, with a depiction of 
pairs oriented in the receding side of the reference particle. On the other hand, we showed 
that the distribution of close particles recovered its expected fore-aft symmetry at low 
concentrations, but that it still remained anisotropic, with a depletion of pairs oriented 
close to the flow direction. We were able to accurately describe the observed anisotropy in 
the pair distribution function by supposing that permanent doublets were completely ab 



sent. W e then showed that the pair distribution function obtained by 



Batchelor fc Green 



1 1972fcr i in the dilute limit, gBcif), accurately follows the simulation results over a wide 
range of r, including the large increase in the probability of finding pairs of spheres near 
contact, corresponding to r 2. However, gsai''') does not take into account the de- 
pletion of permanent doublets, and it is therefore unable to capture the behavior of the 
distribution in the limit r ^ 2. In fact, in contrast to the divergent behavior of 5bg(?') for 
r ^ 2, the numerical results suggest that g{r) ~ for r less than the minimum distance 
of approach between two spheres in the region of open trajectories (rmin 4 x 10~^). 

For the velocity fiuctuations, we showed that, for an isotropic configurational proba- 
bility of particles surrounding a reference sphere located at r, P{CN\r), the temperature 
tensor is diagonal and that the temperatures in the plane of shear are equal. Moreover, 
we showed that in the dilute limit, the temperature components are proportional to the 
volume fraction, and that the temperature in the plane of shear is larger than that in 
the vorticity direction. Then, by averaging the velocity fluctuations originated in the 
hydrodynamic interactions between two spheres, weighted by 5bg(''): Sind neglecting to 
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the first approximation the effects of the permanent doublets, we computed the tem- 
perature tensor in the dilute limit, and found good agreement with the results of our 
numerical simulations, even for moderately concentrated suspensions. Furthermore, we 
were able to accurately reproduce the whole velocity autocorrelation function in both 
transverse directions, on the basis of only two-particle hydrodynamic interactions. In 
contrast, larger discrepancies were found between the corresponding results for the an- 
gular velocity fluctuations and those obtained in the numerical simulations. However, in 
this case we provided a rough estimate for the lower bound of the fluctuations in the 
dilute limit. 

In order to further investigate the effects of the microstructure on the temperature 
tensor, we performed numerical computations in which we calculated the velocity fluctu- 
ations for a hard-sphere distribution of particles subject to the same simple shear flow. 
We also calculated the asymptotic behavior in the dilute limit, using a uniform pair dis- 
tribution function 5hs('') = 1, and obtained an excellent agreement with the numerical 
results for all the linear and angular velocity fluctuations. 

In addition to the temperature tensor, we presented the full probability distribution of 
the velocity fluctuations for both the linear and the angular velocities, in all directions 
and for three different volume fractions. We observed different functional forms as the 
concentration decreases, from a Gaussian to an Exponential and finally to a Stretched 
Exponential form. 

Finally, we presented a simple correction term, which only depends on the mean angular 
velocity of the spheres in the vorticity direction, which enhances the interpretation of the 
LDV measurements at intermediate and low volume fractions. 
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